//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testLegendrePoly.cc
 * \author Jeremy Roberts
 * \date   11/05/2010
 * \brief  test of class LegendrePoly
 * \note   Copyright (C) 2010 Jeremy Roberts
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date:: 2010-11-08 20:29:53 -0500 (Mon, 08 Nov 2010) $:Date of last commit
//---------------------------------------------------------------------------//

#include <cmath>
#include <iostream>
#include "petscvec.h"
#include "petscsys.h"
#include "petscksp.h"
#include "../diff2d/LegendrePoly.hh"
#include "../linalg/typedefs.hh"
using namespace std;
//---------------------------------------------------------------------------//

int main(int argc, char *args[])
{
    // initialize petsc
    PetscInitialize(&argc,&args,PETSC_NULL,PETSC_NULL);
    cout << " TEST CASE FOR LEGENDRE POLYNOMIAL CLASS " << endl;
    cout << " Building two 'partial current' vectors and expanding " << endl;
    cout << " them to a given order in discrete LP's.  " << endl;
	// "test case"
    // -- build the legendre polynomials for a rectangle with:
    // -- ny = 10 vertical fine meshes
    // -- nx = 20 horizontal fine meshe
    // -- fourth order (--> 0, 1, .. 5 )
	integer nx      = 20;  cout << "     nx = " << nx << endl;
    integer ny      = 10;  cout << "     ny = " << ny << endl;
    integer order   = 4;   cout << "  order = " << order << endl;

    LegendrePoly LP;
    LP.buildMe( nx, ny, order );

    // do some test expansions; do a cosine like shape for X and
    // and a gaussian like thing for Y
    Vec PcurX, PcurY;
    VecCreate( PETSC_COMM_SELF, &PcurX );
    VecCreate( PETSC_COMM_SELF, &PcurY );
    VecSetSizes( PcurX, PETSC_DECIDE, nx );
    VecSetSizes( PcurY, PETSC_DECIDE, ny );
    VecSetFromOptions( PcurX );
    VecSetFromOptions( PcurY );
    Vec PcurXappx, PcurYappx;
    VecDuplicate( PcurX, &PcurXappx );
    VecDuplicate( PcurY, &PcurYappx );
    scalar pcX[nx], pcY[ny];
    integer idxx[nx], idxy[ny];

    scalar dx = 2.0/nx;
    for ( integer i = 0; i < nx; i++ )
    {
        idxx[i] = i;
        pcX[i]  = cos( 3.14159265*dx*(0.5 + i)/2.0 );
    }
    cout << " par. current along x: cos(pi*dx*(0.5+i)/2.0), dx=2.0/nx " << endl;
    scalar dy = 1.0/ny;
    for ( integer i = 0; i < ny; i++ )
    {
        idxy[i] = i;
        pcY[i]  = sqrt( dy*(0.5 + i)/1.0 );
    }
    cout << " par. current along y:  sqrt( dy*(0.5+i)/1.0, dy=1.0/ny " << endl;

    // build the vectors
    VecSetValues( PcurX, nx, idxx, pcX, INSERT_VALUES );
    VecAssemblyBegin( PcurX );  
    VecAssemblyEnd( PcurX );
    VecSetValues( PcurY, ny, idxy, pcY, INSERT_VALUES );
    VecAssemblyBegin( PcurY ); 
    VecAssemblyEnd( PcurY );

    // get the expansion coefficients
    scalar coefsX[order+1];
    LP.expandCur( PcurX, coefsX );
    scalar coefsY[order+1];
    LP.expandCur( PcurY, coefsY );

    scalar norm, none = -1.0, zero = 0.0;
    // x
    VecSet( PcurXappx, zero );
    for (integer i = 0; i <= order; i++)
        VecAXPY( PcurXappx, coefsX[i], LP.Px[i] ); // build approximation
    VecAXPY( PcurXappx, none, PcurX );
    VecNorm( PcurXappx, NORM_2, &norm );
    cout << " L2 norm for x approximation  = " << norm << endl;
    // y
    VecSet( PcurYappx, zero );
    for (integer i = 0; i <= order; i++)
        VecAXPY( PcurYappx, coefsY[i], LP.Py[i] ); // build approximation
    VecAXPY( PcurYappx, none, PcurY );
    VecNorm( PcurYappx, NORM_2, &norm );
    cout << " L2 norm for y approximation  = " << norm << endl;


	PetscFinalize();

	cout << "  ... all finished. now check for errors! " << endl;

    return 0;
}

// example, NN = 9, order = 7
// 1.00000  1.00000  1.00000  1.00000  1.00000  1.00000  1.00000  1.00000
// 1.00000  0.77778  0.33333 -0.33333 -1.22222 -2.33333 -3.66667 -5.22222
// 1.00000  0.55556 -0.16667 -0.83333 -0.94444  0.16667  3.33333  9.55556
// 1.00000  0.33333 -0.50000 -0.73810  0.16667  1.83333  2.00000 -4.66667
// 1.00000  0.11111 -0.66667 -0.28571  1.00000  1.00000 -2.66667 -6.22222
// 1.00000 -0.11111 -0.66667  0.28571  1.00000 -1.00000 -2.66667  6.22222
// 1.00000 -0.33333 -0.50000  0.73810  0.16667 -1.83333  2.00000  4.66667
// 1.00000 -0.55556 -0.16667  0.83333 -0.94444 -0.16667  3.33333 -9.55556
// 1.00000 -0.77778  0.33333  0.33333 -1.22222  2.33333 -3.66667  5.22222
// 1.00000 -1.00000  1.00000 -1.00000  1.00000 -1.00000  1.00000 -1.00000


//---------------------------------------------------------------------------//
//                 end of testLegendrePoly.cc
//---------------------------------------------------------------------------//

